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00 ■ ABSTRACT 

CN ; Feedback processes by active galactic nuclei (AGN) appear to be a key for under- 

standing the nature of the very X-ray luminous cool cores found in many clusters 
' r~ I ' of galaxies. We investigate a numerical model for AGN feedback where for the first 

pi |. time a relativistic particle population in AGN-inflated bubbles is followed within a 

full cosmological context. In our high-resolution simulations of galaxy cluster forma- 
tion, we assume that BH accretion is accompanied by energy feedback that occurs in 
two different modes, depending on the accretion rate itself. At high accretion rates, a 
^ ' small fraction of the radiated energy is coupled thermally to the gas surrounding the 

quasar, while in a low accretion state, mechanically efhcient feedback in the form of 
hot, buoyant bubbles that are inflated by radio activity is considered. Unlike in previ- 

(N ■ 

ous work, we inject a non-thermal particle population of relativistic protons into the 
AGN bubbles, instead of adopting a purely thermal heating. We then follow the sub- 
sequent evolution of the cosmic ray (CR) plasma inside the bubbles, considering both 
its hydrodynamical interactions and dissipation processes relevant for the CR popula- 
. tion. This permits us to analyze how CR bubbles impact the surrounding intracluster 

medium, and in particular, how this contrasts with the purely thermal case. Due to the 
different buoyancy of relativistic plasma and the comparatively long CR dissipation 
timescale we find substantial changes in the evolution of clusters as a result of CR 
feedback. In particular, the non-thermal population can provide significant pressure 
support in central cluster regions at low thermal temperatures, providing a natural 
explanation for the decreasing temperature profiles found in cool core clusters. At the 
k>( \ same time, the morphologies of the bubbles and of the induced X-ray cavities show 

, a striking similarity to observational findings. AGN feedback with CRs also proves 

• efficient in regulating cluster cooling flows so that the total baryon fraction in stars 

becomes limited to realistic values of the order of ~ 10%, more than a factor of 3 re- 
duction compared with cosmological simulations that only consider radiative cooling 
and supernova feedback. We find that the partial CR support of the intracluster gas 
also affects the expected signal of the thermal Sunyaev-Zel'dovich effect, with typical 
modifications of the integrated Compton-y parameter within the virial radius of the 
order of - 10%. 

Key words: methods: numerical - black hole physics - cosmic rays - galaxies: 
clusters: general - cosmology: theory 
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1 INTRODUCTION 

There is a growing body of observational evidence indicat- 
ing that many galaxy clusters harbor in their centre super- 
massive black holes (BHs), which interact in a complex way 

* E-mail: deboras@mpa-garching.mpg.de 



with their surroundings. Recently, spectacular images from 
the XMM-Newton and Chandra X-ray telescopes (e.g. Mc- 
Namara et al., 2005; Fabian et al., 2006; Forman et al., 2007) 
have shed some light on this intricate interaction, showing 
clear imprints of central active galactic nuclei (AGN) in the 
intracluster medium of their host clusters. In particular, 
by comparing observations obtained at radio wavelengths 
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with X-ray images, it is often seen that there are X-ray 
depressions in the innermost cluster regions which spatially 
correspond to significant radio emission (e.g. Owen et al., 
2000; Blanton et al., 2001; Clarke et al., 2005). The most 
plausible candidate for inducing these features is the cen- 
tral BH, which is generating jet-inflated radio lobes, often 
simply dubbed 'bubbles'. 

There has been considerable theoretical effort to under- 
stand the interplay of AGN-driven bubbles with the intr- 
aclustcr medium (ICM) (e.g. Churazov ot al., 2001; Quilis 
et al., 2001; Briiggen & Kaiser, 2002; Ruszkowski & Begel- 
man, 2002; Dalla Vecchia et al., 2004; Omma et al., 2004; Si- 
jacki & Springel, 2006a,b; Sijacki et al., 2007). Most of these 
studies have, for simplicity, considered bubbles that are filled 
with hot thermal gas, and by following their evolution with 
time, they tried to constrain how much heat AGN bubbles 
can deliver to the surrounding medium. However, the syn- 
chrotron and inverse Compton emission detected from the 
radio lobes suggests that they contain a significant amount 
of relativistic electrons and that they arc permeated with 
magnetic fields. Besides relativistic electrons and magnetic 
fields, it is also very plausible that an important part of the 
pressure support in the bubbles stems from relativistic pro- 
tons, especially if the AGN jets are heavy. To date it is how- 
ever not clear, both from am observational and a theoretical 
point of view, what the precise composition of the bubbles, 
the relative mixture of hot gas, and non-thermal particle 
components really is. Also, we have no detailed knowledge 
about the strength and configuration of the magnetic fields 
within the bubbles. 

Regarding the contribution from hot thermal gas, the 
work of Mazzotta et al. (2002) suggested that the bubble 
in the MKW 3s galaxy cluster may have a predominantly 
thermal origin, while Schmidt et al. (2002) found that in 
the Perseus cluster a thermal gas component cooler than 
~ 11 keV filling the bubbles seems ruled out. Moreover, it is 
not clear whether entrainment of thermal gas by the radio jet 
is effective and how this should depend on the BH properties. 
A study of a larger sample of galaxy clusters with central 
X-ray cavities by Dunn et al. (2005) hints towards a rather 
complex picture, leaving several possibilities still open, both 
for significant pressure support by non-thermal protons or 
by hot thermal plasma, and for scenarios where the magnetic 
fields axe filamentary at least in some of the cases. 

Bearing in mind the outlined uncertainties it is still 
highly interesting to investigate the possible effects of cosmic 
rays (CRs) in AGN-inflated bubbles on cool core clusters, 
as considered already by several authors (e.g. Boehringcr 
& Morfill, 1988; Loewenstein et al., 1991; Mathews & 
Brighcnti, 2007; Guo & Oh, 2008; Sanders & Fabian, 2007; 
Ruszkowski et al., 2008). These studies have highlighted dif- 
ferent channels of CR interaction with the thermal intraclus- 
ter gas and have demonstrated that CRs can be dynamically 
important in galaxy clusters, contributing up to 50% of the 
central cluster pressure (but see Churazov et al. (2007) for 
somewhat lower estimates of CR pressure support of order 
of 10 — 20%). Nevertheless, due to the complicated nature 
of CR physics, most of these works were based on analytical 
or ID numerical calculations, and neglected any cosmolog- 
ical evolution of the host clusters, and hence any possible 
dependence on its dynamical state. On the other hand, cos- 
mological simulations of CRs produced at structure format 



tion shocks by, for example, Miniati et al. (2001); Miniati 
(2002); Ryu et al. (2003); Ryu & Kang (2004) have high- 
lighted the importance of a realistic cosmological setting for 
a more fateful generation, distribution and following evolu- 
tion of CRs. 

Recently, Enfilin et al. (2007) have proposed a simpli- 
fied formalism for the treatment of CR protons that it suit- 
able for implementation and use in self-consistent cosmolog- 
ical codes, as subsequently demonstrated by Jubelgas et al. 
(2008) and Pfrommcr et al. (2006). In these numerical stud- 
ies, the CR source processes considered were restricted to su- 
pernovae (SNe) and structure formation shock waves. While 
it was found that these CR sources can affect the interstel- 
lar medium of low mass galaxies significantly, they did not 
prevent the overcooling in the centres of galaxy groups and 
clusters. However, the possible impact of CR protons gener- 
ated by central AGN has not been explored in these works, 
and this forms the primary objective of our current study. 

This paper is organized as follows. In Section 2 we out- 
line the methodology wc have adopted to simulate CR bub- 
bles and to follow their cosmological evolution. In Section 3, 
we present test runs performed for isolated halo simulations 
in order to analyze our numerical model in detail, while the 
bulk of our results from cosmological simulations of galaxy 
cluster formation is described in Section 4. Finally, we dis- 
cuss our findings and draw our conclusions in Section 5. 



2 METHODOLOGY 

In this study we use an improved version of the mas- 
sively parallel TreePM-SPH code GADGET-2 (Springel, 
2005; Springel et al., 2001). The SPH formulation adopted 
in the code manifestly conserves both energy and entropy 
even in the presence of fully adaptive smoothing lengths, 
as implemeiitcd by Springel & Hernquist (2002). In addi- 
tion to the gravity of dark matter and baryons, and to or- 
dinary hydrodynamics, the code tracks radiative cooling of 
an optically thin plasma of hydrogen and helium, immersed 
in a time-varying, spatially uniform UV background (as in 
Katz et al., 1996). A subresolution multi-phase model for 
the ISM is used to follows star formation and supernovae 
feedback processes (Springel & Hernquist, 2003). Further- 
more, we adopt the model for BH seeding and growth that 
has been suggested by Di Matteo et al. (2005) and Springel 
et al. (2005), modified however by the introduction of a sec- 
ond mode of AGN feedback to model radio activity, as de- 
scribed by Sijacki et al. (2007). These two modes are thermal 
heating from quasars at high BH accretion rates (BHARs), 
and mechanical feedback in the form of hot, buoyant bubbles 
occurring at low accretion rates. 

In this study, we combine the AGN feedback prescrip- 
tion at low BHARs with a model for CR treatment that 
has been developed, implemented, and discussed by EnBlin 
et al. (2007), Jubelgas et al. (2008), and Pfrommcr ct al. 
(2006). This extension of our model allows us to account 
for a non-thermal component permeating the radio lobes. 
In this section, we outline the most important features of 
our combined BH and CR models, and we describe in detail 
the apjjroximations we have adopted to model non-thermal 
particle populations in AGN-driven bubbles. 



© 0000 RAS, MNRAS 000, 000-000 



Simulations of cosmic ray feedback by AGN in galaxy clusters 3 



2.1 Black hole model 

We represent BHs in the code by coUisionless sink particles 
of initially very small mass, and wc allow them to grow via 
gas accretion and through mergers with other BHs that hap- 
pen to get sufficiently close. During the growth of structure, 
we seed every new dark matter halo above a certain mass 
threshold (e.g. 5 x 10^^° K'^Mq) with a central BH of mass 
10''' /i^^Mq, provided the halo does not contain any BH yet. 
The seeding is accomplished on-the-fly by frequently invok- 
ing a parallel friends-of-friends algorithm. Once seeded, BHs 
can then grow by local gas accretion, with an accretion rate 
estimated with the Bondi-Hoyle-Lyttleton formula (Hoyle & 
Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952). We im- 
pose a limit equal to the Eddington rate on the maximum 
allowed BHAR. 

As for the BH feedback processes, we assume that they 
are composed of two physically distinct modes, depending 
on the BHAR itself. At high accretion rates, i.e. above 10~^ 
in Eddington units, we assume that the BH is in a radia- 
tively efficient phase, where a small fraction of the bolomet- 
ric accretion luminosity is coupled thermally to the local 
gas particles around the BH, with an efficiency of 5% and a 
spherical injection kernel. Instead, at low BHARs, we con- 
jecture that the BH growth is characterized by a radiatively 
inefficient accretion flow, where most of the feedback is in 
a mechanical form, manifesting itself by hot buoyant bub- 
bles that rise through the intragroup or intracluster medium. 
We relate the bubble energy content, radius and duty cycle 
with the BH physics, and adopt a somewhat higher effi- 
ciency of mechanical feedback of 20%. We note that Allen 
et al. (2006) recently showed that there exists a tight cor- 
relation between the Bondi accretion rates calculated based 
on observed gas temperature and gas density profiles and 
estimated BH masses, and the actual power emerging from 
these systems in relativistic jets. This lends observational 
support to our simple estimates of the BHARs. 



2.2 CR implementation 

We represent the CR population in each gaseous fluid ele- 
ment by a relativistic population of protons, which we ap- 
proximately describe with an isotropic power-law distribu- 
tion function in momentum space. In the simple formalism 
adopted here, this distribution function is fully deflned by 
the power law slope a, its normalization C, and a dimcn- 
sionless low momentum cut-off q: An/dp = Cp^"6{p — q). 
Here 9{x) denotes Heaviside step function and the dimen- 
sionless momenta are expressed in units of nip c. We assume 
for clarity that the only source of CRs is given by the ra- 
diatively inefficient accretion mode of BHs, while we neglect 
other contributions from supernovae or particle acceleration 
at cosmic structure formation shocks, which have been ex- 
tensively discussed in previous work (Jubelgas et al., 2008; 
Pfrommer et al., 2006, 2007, 2008; Pfrommcr, 2008). We in- 
clude CR loss processes in the form of thcrinalization by 
Coulomb interactions and losses due to hadronic interac- 
tions. We expect that within dilute bubbles these loss mech- 
anisms should play a subdominant role compared to the 
pressure loss due to adiabatic expansion of the buoyantly 
rising bubbles. Therefore, wc believe that our simulations 
also are representative, at legist qualitatively, of the case of 



a radio plasma that is dominated by a relativistic electron- 
positron population. 

Even though in the present implementation of CR pro- 
cesses in GADGET-2, accounting for CR diffusion is in prin- 
ciple possible, in this study we refrain from trying to model 
it because this treatment of diffusion is only isotropic at 
present. Nevertheless, during the rise of a bubble, CR dif- 
fusion should be significantly suppressed perpendicular to 
the magnetic field lines that drape around the bubbles, only 
allowing CRs to diffuse efficiently in the wake of the bub- 
bles (sec e.g. Sanders & Fabian, 2007; Ruszkowski et al., 
2008). Thus, isotropic CR diffusion would be a rather poor 
representation of this process, and neglecting diffusion alto- 
gx^tlier probably mimics the possible influence of magnetic 
draping effects with higher realism. Note, however, that the 
impact of magnetic fields on CR diffusion out of the bubbles 
is still poorly understood and is a matter of debate, hence a 
detailed treatment of this issue is beyond the scope of this 
work. 

2.3 A model for CR bubbles 

In the numerical framework we adopted for the description 
of CRs, we need to determine the power law slope a of the 
momentum spectrum representative for our system, its injec- 
tion cut-off ginj , and the fraction /cr of the energy released 
by the BH that actually goes into the CRs. Once these val- 
ues have been chosen we can follow the bubble's evolution 
and the loss processes of the CR component self-consistently, 
and compare with the case where the AGN-driven bubbles 
are purely thermal. 

For the slope a we have tested values from 2.1 to 2.4, 
where a steeper power-law slope corresponds to a distribu- 
tion with more low-energy CRs. These can thermalize faster 
and have a less relativistic equation of state that is closer 
to the 'harder' thermal case. We note that this range of val- 
ues for the spectral slope is consistent with observational 
findings for the electron population in FRI sources, espe- 
cially considering young systems (e.g. Birzan et al., 2004; 
Dunn et al., 2005). Similarly to the case in which CRs are 
produced by supernovae (.lubclgas et al., 2008) we estab- 
lish the injection cut-off ginj below which the CR energy is 
instantly thermalized. Equating the injection and loss time 
scale, we solve for the injection cut-off ginj using equations 
(18) and (19) from Jubelgas et al. (2008) and by assuming 
an initial value of the intrinsic injection cut-off q-mit ~ 1 for 
simplicity. Our results do not depend on a particular choice 
for ginit, provided ginit < 1 since Coulomb losses rapidly re- 
move the low energy part of the CR spectrum, which gets 
almost instantly thermalized. 

Finally, at a fixed mechanical feedback efficiency, we can 
choose whether to fill the AGN-inflated bubbles with rela- 
tivistic gas exclusively, or to allow some of the hot thermal 
gas to permeate them as well, by regulating the parameter 
for the energy fraction going into CRs, /cr. Note that, if 
Coulomb losses are important some fraction of initial energy 
delivered by AGN will be thermalized instantly as described 
above, such that /cr can regulate only the remaining frac- 
tion of energy that in principle is available for injection into 
the CR population directly. So far there is no clear consensus 
in the observational literature (e.g. Mazzotta et al., 2002; 
Schmidt et al., 2002; Sanders & Fabian, 2007; Simionescu 
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et al., 2008) whether there is a significant thermal compo- 
nent inside radio bubbles or not, thus one of our main aims 
is to understand to what extent the numerical simulations 
predict noticeable changes in our results for different values 
of /cR , which can provide interesting insights with respect 
to this question. 



3 CR BUBBLES IN ISOLATED HALO 
SIMULATIONS 

In this section wc focus on idealized simulations of isolated 
galaxy clusters that harbour an AGN at their centre. The 
initial conditions for the runs consist of a gaseous cluster at- 
mosphere in hydrodynamic equilibrium within a static NFW 
(Navarro et al., 1996, 1997) dark matter halo, as described 
in more detail in Sijacki & Springel (2006a). These systems 
are then evolved under the influence of radiative cooling, star 
formation, BH growth and different types of AGN feedback. 
The simplicity of these simulations relative to a full cos- 
mological setting permits us to unambiguously identify the 
specific signatures of AGN-inflated bubbles that are filled 
with CRs, and to understand the differences that occur with 
respect to purely thermal bubbles. 

3.1 Thermal versus CR content of bubbles 

We start our analysis by considering a 10^^h~^MQ halo 
where in regular time intervals of 10® yr a pair of spher- 
ical bubbles is injected into the innermost cluster region. 
We have fixed the energy content of the bubbles to be 
2.5 X lO^^erg each, and the bubble radius to be 30/i~^kpc. 
We have performed two identical runs, with the only differ- 
ence being that in the first run the energy was injected as 
thermal energy into the bubbles, while in the second run the 
bubbles were filled with CRs instead. 

In Figure 1, the panels from top to bottom show a time 
evolution sequence of this AGN-heated cluster. The pan- 
els on the left-hand side are projected temperature maps of 
the central cluster region for the case with thermal bubbles, 
while the panels of the middle column show the correspond- 
ing region in the run where the bubbles are filled with rel- 
ativistic particles. To more clearly show the location of the 
CR bubbles, the panels on the right-hand side give the ratio 
of the projected CR pressure to the projected thermal pres- 
sure, so that regions that are permeated with a significant 
CR component can be clearly seen. 

It is evident from this figure that the CRs modify both 
the bubble morphologies as well as their dynamics. The most 
striking difference in the temperature maps is that the CR 
bubbles actually appear as X ray cavities, as has been ob- 
served in a number of clusters (e.g. Owen et al., 2000; Blan- 
ton et al., 2001; Nulsen et al., 2005; McNamara et al., 2005; 
Forman et al., 2007; Fabian et al., 2006). This is due to 
the fact that once the bubbles expand rapidly, their pres- 
sure drops more slowly than in the thermal case due to their 
softer equation of state, so that they become larger and ther- 
mally cooler*, as can be seen by comparing in the maps the 

* Note tfiat prior to the CR injection into the bubble particles 
these have already by construction thermal energy content corn- 



characteristic sizes of bubbles in the cases with thermal and 
CR injection. Thus, even though the bubble radius was the 
same at the moment of injection in both cases, the sub- 
sequent evolution leads to a marked difference in the size 
evolution of the bubbles. 

Furthermore, the CR bubbles remain coherent for a sig- 
nificantly longer time than the thermal bubbles, and they 
reach larger cluster-centric distances as well. Support for 
this finding can be seen in the middle panels of Fig. 1, 
where the thermal bubbles that wore still clearly defined 
in the top panel are already shredded when reaching a dis- 
tance of ~ 150/i~^kpc, while the CR bubbles have risen to 
~ 200 /i~^kpc and are still intact at the same time. One rea- 
son for this behaviour lies in the significantly longer cooling 
time of the rclativistic gas inside the bubbles compared with 
the thermal cooling time. Thus, while we injected initially 
the same energy in both runs, the CR bubbles preserve their 
entropy content for a longer time. Another reason lies in the 
enhanced buoyancy of the CR bubbles thanks to their larger 
size. 



3.2 Self-regulated CR bubble heating 

Now we turn to a more detailed analysis that involves a 
growing BH. To this end we introduce at the beginning of the 
simulations a seed BH of small mass at the centre in order to 
see how feedback by CR bubbles fares in establishing self- 
regulated AGN activity. These test runs are analogous to 
the ones presented in Sijacki et al. (2007), where a more 
detailed description of the set-up can be found. We also test 
variations of the CR power law slope from a value of 2.1 
to 2.4, as illustrated in Figure 2, but we fix ginit = 1 and 
/cR = 1. 

First, wc analyze the growth of the central BH as a 
function of time in the top panels of Figure 2. It can be seen 
that the CR bubbles lead to a self-regulated growth of the 
central BH similar to the thermal case, but at the end of the 
simulated time-span, the BH reaches a slightly smaller mass 
in the CR case. However, the detailed time evolution in the 
case of CR bubbles is physically more complex than in the 
thermal case. Initially, for t < O.ltnubbie, the CR bubbles 
are less effective in heating the ICM, and consequently the 
central BH is growing more efficiently. During these early 
stages of BH growth, the average CR pressure in the clus- 
ter central region is smaller than the central pressure in the 
runs with thermal bubbles. Hence, during this early phase 
the combined effect of Coulomb and hadronic losses as well 
as the higher compressibility of composite CR plus ther- 
mal gas leads to a less efficient CR bubble feedback, as will 
be discussed in more detail in Section 3.3. However, even 
though this allows the gas to cool somewhat more efficiently 
towards the centre the amount of stars formed is roughly 
the same in the thermal and relativistic cases. 

For t > O.ltHubbic, the star formation rate begins to 
be more suppressed in the run with CR bubbles, because at 
this stage Pgr starts to be comparable to the local Pth. In 

parable to the surrounding ICM particles to which wc add non- 
thermal energy released by the AGN. This thermal component 
inside the bubbles is in most of the cases negligible with respect 
to the CR component. 
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Figure 1. Projected mass- weighted temperature and pressure maps of tlie central cluster region at three different evolutionary stages, at 
times: O.OTtjjubblo (top panels), 0.12 tjj^jji,!;, (middle panels), and 0.24 ^JJ^^,^,l(, (bottom panels). The left-hand panels show temperature 
maps in the run where thermal bubbles have been injected in regular time intervals. The central panels illustrate a simulation where the 
bubbles were initially filled with relativistic particles instead. For this run, the regions where the CR pressure is significant with respect 
to the thermal pressure are shown in the right-hand panels. It can be clearly seen that the bubble dynamics, coherence and maximum 
cluster-centric distance reached are affected by the presence of a relativistic component filling the bubbles. 



fact, for the run with a = 2.1, this transition occurs some- 
what earlier than in the simulation with a steeper a, given 
that for Q = 2.1 the cosmic ray pressure reaches significant 
values sooner. It can be seen that for t > O.ltHubbic the 
BHAR starts to oscillate dramatically as a consequence of 
the presence of a dominant relativistic particle component. 
A similar behaviour has been observed in the case of star 
formation in small dwarf galaxies in the work by Jubelgas 
et al. (2008). When a bubble is injected into the ICM and 
its Pgr is comparable or higher than the local Pth it will 
provide an extra pressure support to the ICM. While rising 



buoyantly towards the cluster outskirts and heating the sur- 
rounding gas, this dramatically reduces the amount of gas 
that can be accreted by the central BH. Successively, CR 
pressure will get dissipated and gas will start to cool again 
leading to higher BHAR values, until another bubble is trig- 
gered, and the cycle is repeated. Note that at late times due 
to this behaviour of CR bubbles there is some residual star 
formation activity and the BH growth is not as efficiently 
terminated as it is the case for thermal bubbles. 

The bottom panels of Figure 2 show how the profiles 
of gas density, temperature and entropy are affected by CR 
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Figure 2. The top panels show the BH mass growth, BHAR in Eddington units, and star formation rate (SFR) as a function of time 
for four different runs: black lines denote the case with cooling and star formation only, blue lines show results where the AGN-driven 
bubbles are thermal, while the green and red lines refer to the case of CR bubbles. The run with a steeper spectral index of a = 2.4 is 
represented with green lines, while the case of a = 2.1 is drawn with red lines. The bottom panels show radial profiles of gas density, 
mass-weighted temperature, and entropy, for the same set of runs and with the same colour-coding. The vertical dotted lines denote the 
gravitational softening length and the virial radius of this halo. 



bubbles. The radial temperature profile is always decreasing 
towards the central regions in the runs with CR bubbles, 
reflecting the fact that a significant fraction of non-thermal 
pressure builds up there. However, the gas density can ei- 
ther be increased or reduced relative to the thermal case, 
depending on whether Pgr is higher than Pth at a given 
instant. Finally, the gas entropy is also somewhat reduced 
in the innermost regions, mainly being driven by the varia- 
tions in temperature. Overall, the simulation with a flatter 
CR spectrum shows a more significant imprint of CRs on the 
ICM, as expected due to the more pronounced high-energy 
tail of the CR momentum distribution, so that a significant 
level of Pgr can be maintained for a longer time interval. 



3.3 Importance of Coulomb and hadronic losses 

There are two primary dissipation processes for the CR par- 
ticle population. One is given by individual electron scat- 
terings in the Coulomb fields of CR particles and by small 
momentum transfers through excitations of plasma oscilla- 
tions, which is especially efficient in the low momentum tail 
of the CR population (Gould, 1972). The other important 
CR energy loss process is due to inelastic collisions of CRs 



on nucleons, where hadronic interactions produce pious that 
are then ultimately dissipated into radiation as well as ther- 
mal heat (e.g. Enfilin et al., 2007). Additionally, the pro- 
cess of MHD wave-mediated cosmic-ray heating (Loewen- 
stein et al., 1991) could be relevant for the CR particle pop- 
ulation within the bubbles (see e.g. Guo & Oh, 2008), that 
we do not consider in this work. 

In order to gauge the relative importance of these loss 
processes we have performed two additional test runs where 
we have switched off CR losses altogether, or alternatively 
only allowed for Coulomb cooling. Based on these runs we 
find that in the case where bubbles are filled exclusively with 
CRs and where no diffusion is allowed, hadronic losses are 
not very important during the whole simulated timespan. 

Instead, we find that Coulomb cooling is especially im- 
portant in the early phases of BH growth. Here the Coulomb 
cooling time of the injected CR particles turns out to be 
rather short such that the CR energy is quickly thermalized 
into the ICM. Thus, in the run where only Coulomb cool- 
ing is allowed the BH is accreting mass as efficiently in the 
beginning as in the case with thermal feedback. However, 
for t > 0.1 f Hubble, the CR pressure inside of the bubbles 
builds up and the characteristic time of Coulomb losses in- 
creases. As a result, the CRs thermalize very little of their 
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Figure 3. Maps of the central region of the 3676 galaxy cluster at 2 = 0.1. In the upper row, a mass- weighted temperature map is shown 
in the left-hand panel, while the right-hand panel illustrates the mass-weighted spatial distribution of PcR relative to the total pressure. 
In the lower row. X-ray luminosity maps are illustrated. The right-hand panel gives an unsharp masked luminosity map, allowing for 
easier bubble identification. 



energy content for t > O.lSfHubbie- Interestingly, in the case 
where no losses were allowed for the non-thermal population 
BH growth is somewhat faster during the very first part of 
the simulation than in the case of thermal bubbles. Here, 
CRs cannot thermalize part of their energy into the sur- 
rounding ICM and heat it, and since the composite of non- 
thermal and thermal pressure is more compressible, AGN 
feedback becomes somewhat less efficient than the purely 
thermal bubbles in the early stages of BH growth. 

We have further explored the case where only one part 
of the available mechanical energy is delivered into CRs, 
while the other remained in thermal form. For values of /ch, 
as low as 0.1, the influence of CRs is barely distinguishable 



from the purely thermal case. However, if we select /cr of 
order of 0.5 or higher the composite of CRs plus thermal gas 
in the bubbles shows interesting properties. In particular, 
for t < 0.1 tHubbio, the non-thermal component inside of the 
bubbles only weakly increases the gas compressibility, and 
these composite bubbles essentially behave as pure thermal 
bubbles. On the other hand, for t > O.ltnubbio, the cosmic 
ray pressure Pea inside of the composite bubbles reaches 
high values, so that CR processes become much more rele- 
vant. At this stage, the properties of composite bubbles and 
their feedback effects become very similar to the case of pure 
CR bubbles. 
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4 COSMIC RAY BUBBLES IN SIMULATIONS 
OF CLUSTER FORMATION 

In this section we consider self-consistent cosmological simu- 
lations of galEixy cluster formation subject to AGN heating. 
Prom a cosmological box of size 479 ft.~^Mpc on a side that 
has been evolved adopting a ACDM cosmology (Yoshida 
et al., 2001; Jenkins et al., 2001), a galaxy cluster with a 
mass of ^ lO"'* Ii^^Mq has been selected (the 'g676' clus- 
ter), that is fairly rela^xed at the present epoch. Using the 
Zoomed Initial Condition (ZIC) technique (Tormen et al., 
1997), new initial conditions for this cluster have been con- 
structed by resampling its Lagrangian region at a higher 
numerical resolution (Dolag, 2004). In Sijacki et al. (2007), 
we have presented several simulations of this cluster with 
diflFerent physics, including runs with cooling and star forma- 
tion processes only, as well as simulations that analyzed the 
impact of our AGN feedback model with thermal bubbles. 
Here, we discuss analogous simulations that focus on the 
role of CRs in AGN-inflated bubbles. Further details about 
the numerical setup can be found in Sijacki et al. (2007). 
In our new simulations, wc keep all numerical parameter 
choices (like gravitational softening lengths, number of SPH 
smoothing neighbours, etc.) exactly the same in order to 
clearly pin down the impact of the non-thermal nature of 
the bubbles. 



4.1 Signatures of CR bubbles 

We start our analysis by trying to identify a link between 
the CR bubble properties and specific features they induce 
in the ICM. In Figure 3, we show the central region of our 
most massive cluster at z = 0.1. The thickness of the pro- 
jection slice is 300ft~^kpc, and the different panels give the 
projected mass-weighted temperature map (upper left-hand 
panel), the PcR-to-total pressure ratio (upper right-hand 
panel), the X-ray luminosity map (lower left-hand panel) 
and an unsharp masked luminosity map (lower right-hand 
panel). The unshaxp masked luminosity map has been pro- 
duced by dividing the original luminosity map by a version 
of it that has been smoothed on ~ 180/i~^kpc scale. In the 
upper right-hand panel, the position of CR bubbles can be 
clearly identified. In the very centre of the cluster there is 
significant cosmic ray pressure, corresponding to very re- 
cent BH activity which is producing two small bubbles. 
Moreover, upwards from the central region there is a large 
plume filled with CR protons extending from ~ 200/i^^kpc 
to ^ 500/i~^kpc along the z-axia. Yet another plume is 
clearly visible in the lower region, propagating towards the 
lower-left corner, and having a total length of more than 
~ 200ft^^kpc. The temperature map reveals asymmetries 
along the direction of big CR plumes, and a close corre- 
spondence to the central small bubbles, where the temper- 
ature is somewhat reduced. However, it is in the luminosity 
maps that one can clearly identify CR bubbles, as depres- 
sion regions in X-ray emissivity. The bubble morphologies 
and the reduction in X-ray emissivity in the regions that 
are permeated with CR particles are in very good qualita- 
tive agreement with a number of observational findings (c.g 
Owen et al., 2000; Nulsen et al., 2005; Forman et al., 2007; 
Sanders & Fabian, 2007), suggesting that the presence of a 
non-thermal component inside of the bubbles is essential for 



more faithful reproduction of the observational results and 
their understanding. 

In Figure 4, we analyze the CR bubble morphologies at 
diflFerent cosmic times. To this end we plot the projected 
CR energy per unit mass in the x-z plane within a cu- 
bic box of size 4000ft^^kpc on a side, at three different 
redshifts, z = 1.0, 0.6 and 0.2. The left-hand panels refer 
to the run performed with a = 2.1, while the right-hand 
panels are for a = 2.4. The black dots over-plotted on the 
maps mark the positions of BH particles more massive than 
1.5 X 10^ h~^MQ. Note that we have deliberately chosen not 
to vary the injection axis of the bubbles with time (which 
is always along the 2;-axis). This has been done in order to 
assess to what extent environmental effects like bulk mo- 
tions, local ovcrdcnsity fluctuations, or merger events, can 
displace bubbles from their initial injection position. At the 
first glance it can be noticed that the evolved CR bubble 
morphologies are rather complex, clearly lacking a strong 
axial symmetry, oven though they were injected in a sym- 
metric fashion. This emphasizes the important role played 
by the cosmological environment in modifying the bubble 
shapes. 

The three epochs shown in Fig. 4 have been selected 
such that they correspond to merging events of the central 
BH with a smaller mass BH that once belonged to a satel- 
lite halo. Thus, at these epochs, the BH activity in the cen- 
tral region is rather vigorous, and several bubble episodes 
can be seen, from the central ones that are just injected, to 
the more peripheral ones which have risen and are now de- 
tached, but still have significant pressure support from CR 
protons. Interestingly, a,t z = 1.0, in the vicinity of the most 
meissive cluster in the lower-right part, there is a smaller 
halo containing a central massive BH which during a sim- 
ilar time span has also generated two bubbles. At z = 0.6 
this smaller halo reaches ~ 1/3 of the central cluster mass 
and merges with it violently. The central BHs of these halos 
merge as well, increasing the mass of the most massive BH 
by ~ 10'' H^^Mq. It can be seen that during this merger 
event the relic CR bubbles are significantly perturbed and 
displaced towards the right, in the direction from which the 
smaller halo was falling in. 

Finally, by comparing the left-hand panels to the right- 
hand panels in Fig. 4, it can be seen that CR bubbles with 
a steeper power-law spectrum have in general similar mor- 
phologies to the a = 2.1 case, but they show a somewhat 
smaller size and shorter survival time, as expected given 
that Coulomb losses become more important for a steeper 
spectrum. 

4.2 Impact of CR bubbles on their host haloes 

4-2.1 Radial profiles 

In Figure 5, we show radial profiles of gas density (first row), 
mass- weighted temperature (second row), entropy (third 
row), and thermal pressure (last row) at three different 
epochs, z = 0.5, 0.2 and 0.0. The black continuous lines 
are for the run with cooling and star formation only, blue 
dashed lines denote the run with thermal bubbles (as in Si- 
jacki et al., 2007), while the green dot-dashed (a = 2.4) 
and red triple dot-dashed lines (a — 2.1) show the case 
with CR bubbles. It can be seen that the gas density is 
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Figure 4. Projected mass-weighted CR energy per unit mass in boxes of 4000?t~^kpc on a side at three different redshifts: z = 1.0 
(top), 0.6 (middle) and 0.2 (bottom). The left-hand panels are for the run with a = 2.1, while the right-hand panels show the case of 
a = 2.4. The selected redshifts correspond to highly perturbed states of the cluster when it is undergoing merging events. It can be seen 
how environmental effects influence the morphologies and dynamics of CR bubbles. 
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Figure 5. Radial profiles of gas density (first row), mass- weighted temperature (second row), entropy (third row) and thermal pressure 
(last row, in internal code units) of the g676 galaxy cluster at z = 0.5, 0.2 and 0.0, respectively. Solid black lines illustrate the case 
without AGN heating. The dashed blue lines are for simulations where AGN feedback is included with purely thermal bubbles, as in 
Sijacki et al. (2007). The runs with CR filled bubbles are represented with green dot-dashed lines and red triple dot-dashed lines, for 
spectral indices of a = 2.4 and a = 2.1, respectively. The vertical dotted lines denote the gravitational softening length and the virial 
radius of this cluster. The dash-dotted lines in the second row of panels show the slope of the central temperature profiles of the cool 
core clusters found by Sanderson et al. (2006), while the dash-dotted lines in the third row illustrate the entropy scaling with radius, 
i.e. S oc r^'^. In the last row we also indicate the mean ratio of CR-to-total pressure in each bin, which is shown with green diamonds 
(a = 2.4) and red triangles (a = 2.1). 
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affected by the presence of relativistic particles in the bub- 
bles, being roughly comparable to the simulation where the 
bubbles are purely thermal. However, the gas temperature is 
altered substantially and in general is always reduced in cen- 
tral regions, as wc have also found for our isolated cluster 
simulations. This is simply caused by the fact that signif- 
icant pressure support in the central cluster regions is in 
a non-thermal form, which decreases the thermal pressure 
needed to counteract the gravity exerted by the cluster's 
potential. In Sijacki ct al. (2007) wo have already shown 
that a sclf-rcgulatcd AGN feedback mechanism is required 
in order to obtain realistic temperature profiles that can re- 
produce the decline towards the innermost regions that has 
been observed in a number of cool core clusters. Having a 
significant population of non-thcrrrial particles in the cen- 
tral bubbles further helps in this respect, and provides a 
natural explanation for the observed shape of the temper- 
ature profiles. Thus, the presence of a relativistic particle 
population inside AGN-inflated bubbles is not only essen- 
tial for reproducing the observed morphologies of bubbles 
and X-ray cavities in the hot ICM, but it also improves the 
temperature distribution of the simulated galaxy clusters. 
However, our results also show that the AGN feedback effi- 
ciency and the assumed fraction of the energy in high energy 
CR protons are to a certain extent degenerate, suggesting 
that a detailed comparison with observations is needed to 
constrain both of these parameters individually. 

The CR bubble feedback also steepens the entropy pro- 
file and brings it into closer agreement with X-ray obser- 
vations, especially at intermediate radii from ~ 20/i^^kpc 
to ~ 300 h~^kpc. Further information about the role of CR 
bubbles can be obtained from the last row of panels in Fig- 
ure 5, where we plot radial profiles of the mass-weighted 
thermal pressure for all the runs we performed. The symbols 
give the mean ratio of CR-to-total pressure in each spherical 
bin. It can be seen that the CR pressure is most relevant in 
the central regions, being comparable to the thermal pres- 
sure for r < 50/i^^kpc at all throe epochs considered. At 
z = and inside 100/i~^kpc, the mean CR pressure con- 
tributes about 37% of the total pressure for a spectral slope 
of Q = 2.1 (for the case of a = 2.4 this number is ~ 23%), 
while at the virial radius, the non-thermal pressure support 
drops to 14% (or 10% for a = 2.4). Recall however that in 
this study we have excluded other potential sources of CRs, 
namely supernovae and cosmic structure formation shocks. 
They arc expected to further increase the overall CR pres- 
sure content, and may also affect the relative contributions of 
iion-tlicrmal to thermal pressure as a function of radius (for 
further details on these issues see Pfrommer et al., 2007). In- 
terestingly, given that the CR pressure contribution is com- 
parable to the thermal pressure and that the CR pressure 
gradient is negative, the central cluster gas might be convec- 
tively unstable as proposed by Chandran & Dennis (2006) 
and explored in the case of SNe injected CRs by Pfrommer 
et al. (2007). Actually, since most of the CR pressure is due 
to CRs residing in the buoyantly rising bubbles, one can ar- 
gue that the system indeed exhibits a CR-driven convection. 
An instability did not need to develop, since the CRs in the 
ICM are directly injected into an unstable configuration of 
a bubble. 

The total X-ray luminosity within the virial radius in 
the case of CR bubbles with a = 2.4 lies in between the 



simulations performed without AGN feedback and the run 
with thermal bubbles. This increase of the X-ray emission 
is a result of the substantial Coulomb losses of the CR pop- 
ulation combined with the softer equation state of the gas 
supported by non-thermal pressure. In the run with a — 2.1, 
Coulomb losses arc less important and the total X-ray lumi- 
nosity is reduced from a value of 1.6 x 10*"^ erg/s for the 
case without AGN feedback to ~ 9.2 x 10''^ erg/s. Also, in the 
simulation with a = 2.1, the mean mass- weighted gas tem- 
perature is reduced, from 1.20 keV to 0.97 keV, indicating 
that CR bubbles may significantly influence galaxy cluster 
scaling relations such as the Lx — T relationship, especially 
at the low mass end. We plan to analyze this question in 
more detail in a forthcoming study. Finally, we note that 
the bulk of the BH growth is not affected rrmcli by the pres- 
ence of a non-thermal component in the bubbles: the BH 
growth is moderately reduced for 2; < 2, so that the most 
massive BH in the simulated volume has its mass reduced 
by less than 20% at « = 0. We also note that in our model 
the central star formation rate and the total stellar mass 
within iivir (see below) are not very sensitive to whether 
the bubbles are filled with thermal or non-thermal energy, 
indicating that feedback by CR bubbles is about equally ef- 
ficient in suppressing star formation in central galaxies as 
thermal ones. 



4-2.2 Baryon fraction and AGN-driven bubbles 

In the left-hand panel of Figure 6, we show the gas fraction 
at 2: = within a given radius as a function of distance 
from the cluster centre. The radius has been normalized to 

i?2()() and the gas fraction is in units of the universal baryon 
fraction adopted in our simulation (i.e. Q,b/Q,o = 0.13). The 
black line is for the run without AGN feedback, the blue 
line is for thermal bubbles, while the green (a = 2.4) and 
red (a = 2.1) lines show the case of CR bubbles. As a generic 
feature of AGN feedback, it can be seen that the gas frac- 
tion is significantly reduced in the central cluster region, 
while it is increased in the cluster outskirts. As expected, 
the central gas fraction is most suppressed in the simula- 
tion where the CR bubbles have a shallow momentum spec- 
trum. On the other hand, the case with a = 2.4 leads to a 
somewhat higher gas fraction with respect to the simulation 
with purely thermal bubbles, due to the combined effects of 
stronger Coulomb cooling and higher gas compressibility. 

In the right-hand panel of Figure 6, we show the total 
baryon fraction (solid lines), and the stellar fraction (dashed 
lines) as a function of radius, using the same colour-coding. 
The amount of stars that form in the central regions is signif- 
icantly reduced by the presence of a supermassive BH, but 
it is not affected by the nature of the bubbles, as discussed 
above. The total amount of stars drops from a high value 
in excess of ~ 30% reached in the simulation without AGN 
heating to a much more realistic value of ~ 10%, which is 
consistent with observational findings (Lin et al., 2003). This 
also endows the central galaxy with a much redder colour, as 
observed. Due to the reduction of the central stellar fraction 
the total baryon fraction is lowered in the central regions as 
well. Only at radii comparable to the turnaround radius, at 
~ 3 X -R200, it reaches the universal value in all runs per- 
formed. 
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Figure 6. The left-hand panel shows the profile of the cumulative gas mass fraction as a function of radius at a = 0. The baryon fractions 
have been normalized to the universal value assumed in our simulations, i.e. 0.13. The black solid line is for the simulation with cooling 
and star formation only, while the remaining lines illustrate results for different runs with AGN feedback: thermal bubbles (blue dashed 
line), or bubbles filled with CRs (green dot-dashed for a = 2.4, and red triple dot-dashed for a = 2.1). The right-hand panel shows the 
total baryon mass fraction (solid lines) and the stellar mass fraction (dashed lines) with the same colour-coding. 
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Figure 7. Ratio of the Compton-j/ parameters in our simulations without AGN feedback and in runs with AGN-driven bubbles at 2 = 0. 
The left-hand panel shows the case of thermal bubbles, the central panel that of CR bubbles with a = 2.4, and the right-hand panel 
the case of CR bubbles with a steeper spectral index of a = 2.1. It can be seen that the Compton-j/ parameter is decreased in the 
central regions due to the bubbles in all runs considered. However, in the cluster outskirts the bubbles tend to increase the Compton-j/ 
parameter, such that the area-integrated Comptonization can be both lower or higher than in the runs without AGN feedback. 



4.2.3 Thermal SZ effect 

Given that CRs provide a significant pressure support in 
the central cluster regions it is interesting to examine how 
they affect the amplitude of the thermal Sunyaev-Zel'dovich 
effectt. In Fi gure 7, we plot the ratio of the Compton-?/ pa- 

t Note that we do not need to consider relativistic corrections to 
the pure thermal Sunyaev-Zel'dovich effect, given that we have 
modelled only CR protons. For a comprehensive study of non- 
thermal bubble imprints on the Sunyaev-Zel'dovich effect see e.g. 
Pfrommer et al. (2005). 



rameter between the runs with and without AGN feedback 
at 2 = 0. The left-hand panel illustrates the case where 
bubbles are purely thermal, while the central {a = 2.4) 
and right-hand panels (a — 2.1) refer to the runs with 
CR bubbles. The Compton-j/ parameter is essentially de- 
termined by the line-of-sight integral of the thermal pres- 
sure. Since AGN-driven bubbles modify the thermal pres- 
sure distribution within the cluster, we observe variations of 
the estimated Compton-?/ parameter. More specifically, the 
Compton-j/ parameter is significantly reduced in the cen- 
tral cluster regions for all three runs with AGN feedback. 
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As a side effect of solving the overcooling problem, the self- 
regulated AGN feedback results in a lower central density 
that is in better agreement with the density profiles derived 
from X-ray measurements. In the cluster outskirts, however, 
the thermal pressure is somewhat increased in the runs with 
AGN feedback. This implies that the gas pressure distribu- 
tion is somewhat "puffed up" in the cluster outskirts, which 
compensates for central gas depletion. Given that the ther- 
mal energy of the cluster is roughly constant and indepen- 
dent of the presence or absence of AGN heating, the area 
integrated Compton-y parameter turns out to be quite ro- 
bust and is not changed much. It is primarily a measure of 
the cluster's gravitational potential. 

Still, this balance is not perfect eis can be deduced by 
calculating the total Compton-y parameter Y within -Rvir. 
In the case of thermal bubbles, we find that the integrated 
Compton-i/ parameter is increased by ~ 10%, while in the 
case of CR bubbles with spectral index of 2.4 this number 
is somewhat smaller, being of order 7%. Interestingly, in the 
case of CR bubbles with a flatter spectral slope the total 
Compton-y parameter is instead lowered by ~ 7%, as a re- 
sult of a substantial reduction of the thermal pressure in the 
innermost regions. These results indicate that the integrated 
Comptonization is sensitive at a level of ^ 10% to the na- 
ture of AGN-driven bubbles and to the adopted feedback 
efficiency. This suggests that the presence of AGN feedback 
contributes to the scatter in the scaling relation between the 
cluster mass and Y, and the magnitude of this effect depends 
on the detailed physics of the AGN feedback. Further studies 
are needed to characterize the functional form of the scatter 
and to quantify a possible mass dependence of this effect. It 
is clear that the intended use of y-measurements for high 
precision cosmology with galaxy clusters will require an ac- 
curate understanding of the detailed physical properties of 
AGN-inflated bubbles. 



5 DISCUSSION 

In this study we have investigated a model for self-regulated 
BH feedback with non-thermal relativistic particles in AGN- 
inflated bubbles. To describe the cosmic ray particles, we 
have adopted the formalism developed by Enfiliu ct al. 
(2007) and Jubelgas et al. (2008) for the treatment of CR 
protons, and combined it with prescriptions for BH growth 
and feedback, outlined in Springel et al. (2005) and Sijacki 
et al. (2007). This allowed us to study the influence of CR 
bubbles on their host clusters in self-consistent cosmological 
simulations of galax:y cluster formation, and in particular, 
to compare with the case where the bubbles are filled with 
purely thermal gas. Our methods also allowed us to gain 
some insight on how strongly the bubble morphologies and 
their propagation through the cluster atmosphere are af- 
fected by a cluster's dynamical state and the gaseous bulk 
motions in the ICM. 

We have found that simulations with CR bubbles much 
more faithfully resemble observational findings compared 
with the results for thermal bubbles. This is mostly driven by 
the softer equation of state of gas that is supported by non- 
thermal pressure, and by the longer dissipation timescale 
of the CR component relative to thermal cooling. As a re- 
sult, CR bubbles can rise to the cluster outskirts and even 



leak into the surrounding intergalactic medium. Still, their 
pressure support becomes gradually less important as they 
move away from cluster central regions. We have also found 
that the bubble morphologies in forming clusters are rather 
complex, especially during major merger events, where bub- 
bles can become significantly perturbed from their initially 
spherical shape and arc displaced from their initial injection 
axis. This implies that based on the spatial position of bub- 
bles alone it is not readily possible to infer whether the AGN 
jet has changed its orientation between two successive out- 
bursts or whether bulk motions in the ICM have displaced 
buoyantly rising bubbles. 

Interestingly, we have found in our cosmological simu- 
lations that neither the BHARs nor the SFRs are very sen- 
sitive to the nature of the bubble feedback. However, the 
temperature distribution in the ICM is affected noticeably, 
given that CR bubbles provide non-thermal pressure sup- 
port in central cluster regions. Above all, this changes the 
galaxy cluster temperature profiles which decline towards 
the innermost regions, being in much better agreement with 
the observed temperature profiles of cool core clusters. Cor- 
respondingly, the thermal pressure distribution is altered as 
well by CR bubbles, with a reduction of the thermal pressure 
in the centre and an increase in the outskirts. This redistri- 
bution of the thermal pressure can modify the integrated 
Compton-j/ parameter at a level of ~ 10%, an effect that 
we expect to contribute to the scatter in the relationship 
between cluster mass and Y. 

We have evaluated how many significant bubble out- 
burst the central BH in our cluster is undergoing during the 
time interval from z = 2 to z = 0, when the BH is in the 
low accretion, 'radio-mode' regime. Based on this we obtain 
an average duty cycle of ~ 10** yrs. Wo note however that 
sudden infiows of gas towards the central BH triggered by 
merging activity - especially at z = 0.6, which corresponds 
to the last major merger of the host cluster - can cause much 
more frequent bubble events than the average value during a 
limited period. There is hence considerable variability in the 
duty cycle related to the merger history of the host halo, an 
issue we plan to investigate in more detail in a forthcoming 
study. 

In our simulations, the sizes of the bubbles arc not de- 
termined by calculating the jet physics from first principles, 
given that we cannot resolve the initial stage of jet forma- 
tion and bubble inflation by a jet even in the most advanced 
statc-of-art cosmological simulations. Instead, we rely on 
simple physical parameterizations and observational guid- 
ance to inject the bubbles with plausible sizes. Certainly 
this is an important limitation of our work and may in part 
be responsible for the rather large bubble sizes of up to a 
few hundred kpc that we find in some cases. On the other 
hand, there are two additional factors that may favor rather 
large bubble dimensions in our current modelling. One is 
that we have kept the injection axis constant with time. In 
some cases, the bubble injection time matches the buoyant 
rise time in such a way that fresh CR protons arc injected 
into an already existing bubbles, making them reach a larger 
size. If this is also realized in nature, it is potentially am- 
biguous whether large bubbles are created by particularly 
powerful outbursts or result from several smaller outbursts 
which combine together. 

However, we note that if CR diffusion is efficient (which 
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we have not considered in this study) it would be less prob- 
able that several bubble outbursts mix together into a large 
X-ray cavity. In this case an old bubble would already dis- 
perse before a new one could have time to reach it and com- 
bine with it. CR diflFusion would also reduce the probability 
of finding a bubble at large enough distance from the cluster 
centre (and would also possibly modify the heating rates), 
thus more sophisticated studies and more stringent obser- 
vational constrains are needed to pin-down the relevance of 
CR diffusion out of the bubbles. 

The second factor which may bias the sizes of our sim- 
ulated CR bubbles high is the fact that observationally only 
relativistic electrons have been detected so far in radio lobes, 
while we have conjectured in our simulations that there 
should bo a spatially extended distribution of CR protons as 
well. Many of our CR bubbles, especially the ones that are 
further away from the central AGN, are expected to have 
a rather aged population of CR electrons, which would not 
have detectable radio emission. In our picture these bubbles 
correspond to the so-called ghost cavities which arc identi- 
fied as depressions in X-ray emission without any notable 
radio emission, or only radio emission at lower frequencies. 
Because of the very low contrast in projected surface bright- 
ness maps, many of these ghost cavities are missed both in 
observations and simulations. It will be an important task 
for future work to better quantify this bias, and to find pow- 
erful observational constraints for the relative content of rel- 
ativistic electrons and protons in AGN-inflatcd bubbles. Our 
simulations certainly suggest that CRs may play a crucial 
role in shaping the central ICM properties and in regulating 
AGN activity in clusters of galaxies. 
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